1 /*
2 * Copyright (C) 2011 The Guava Authors
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 package com.google.common.math;
18
19 import static com.google.common.base.Preconditions.checkArgument;
20 import static com.google.common.base.Preconditions.checkNotNull;
21 import static com.google.common.math.MathPreconditions.checkNoOverflow;
22 import static com.google.common.math.MathPreconditions.checkNonNegative;
23 import static com.google.common.math.MathPreconditions.checkPositive;
24 import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
25 import static java.lang.Math.abs;
26 import static java.lang.Math.min;
27 import static java.math.RoundingMode.HALF_EVEN;
28 import static java.math.RoundingMode.HALF_UP;
29
30 import com.google.common.annotations.GwtCompatible;
31 import com.google.common.annotations.GwtIncompatible;
32 import com.google.common.annotations.VisibleForTesting;
33
34 import java.math.BigInteger;
35 import java.math.RoundingMode;
36
37 /**
38 * A class for arithmetic on values of type {@code long}. Where possible, methods are defined and
39 * named analogously to their {@code BigInteger} counterparts.
40 *
41 * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
42 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
43 *
44 * <p>Similar functionality for {@code int} and for {@link BigInteger} can be found in
45 * {@link IntMath} and {@link BigIntegerMath} respectively. For other common operations on
46 * {@code long} values, see {@link com.google.common.primitives.Longs}.
47 *
48 * @author Louis Wasserman
49 * @since 11.0
50 */
51 @GwtCompatible(emulated = true)
52 public final class LongMath {
53 // NOTE: Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
54
55 /**
56 * Returns {@code true} if {@code x} represents a power of two.
57 *
58 * <p>This differs from {@code Long.bitCount(x) == 1}, because
59 * {@code Long.bitCount(Long.MIN_VALUE) == 1}, but {@link Long#MIN_VALUE} is not a power of two.
60 */
61 public static boolean isPowerOfTwo(long x) {
62 return x > 0 & (x & (x - 1)) == 0;
63 }
64
65 /**
66 * Returns 1 if {@code x < y} as unsigned longs, and 0 otherwise. Assumes that x - y fits into a
67 * signed long. The implementation is branch-free, and benchmarks suggest it is measurably
68 * faster than the straightforward ternary expression.
69 */
70 @VisibleForTesting
71 static int lessThanBranchFree(long x, long y) {
72 // Returns the sign bit of x - y.
73 return (int) (~~(x - y) >>> (Long.SIZE - 1));
74 }
75
76 /**
77 * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
78 *
79 * @throws IllegalArgumentException if {@code x <= 0}
80 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
81 * is not a power of two
82 */
83 @SuppressWarnings("fallthrough")
84 // TODO(kevinb): remove after this warning is disabled globally
85 public static int log2(long x, RoundingMode mode) {
86 checkPositive("x", x);
87 switch (mode) {
88 case UNNECESSARY:
89 checkRoundingUnnecessary(isPowerOfTwo(x));
90 // fall through
91 case DOWN:
92 case FLOOR:
93 return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x);
94
95 case UP:
96 case CEILING:
97 return Long.SIZE - Long.numberOfLeadingZeros(x - 1);
98
99 case HALF_DOWN:
100 case HALF_UP:
101 case HALF_EVEN:
102 // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
103 int leadingZeros = Long.numberOfLeadingZeros(x);
104 long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros;
105 // floor(2^(logFloor + 0.5))
106 int logFloor = (Long.SIZE - 1) - leadingZeros;
107 return logFloor + lessThanBranchFree(cmp, x);
108
109 default:
110 throw new AssertionError("impossible");
111 }
112 }
113
114 /** The biggest half power of two that fits into an unsigned long */
115 @VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L;
116
117 /**
118 * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
119 *
120 * @throws IllegalArgumentException if {@code x <= 0}
121 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
122 * is not a power of ten
123 */
124 @GwtIncompatible("TODO")
125 @SuppressWarnings("fallthrough")
126 // TODO(kevinb): remove after this warning is disabled globally
127 public static int log10(long x, RoundingMode mode) {
128 checkPositive("x", x);
129 int logFloor = log10Floor(x);
130 long floorPow = powersOf10[logFloor];
131 switch (mode) {
132 case UNNECESSARY:
133 checkRoundingUnnecessary(x == floorPow);
134 // fall through
135 case FLOOR:
136 case DOWN:
137 return logFloor;
138 case CEILING:
139 case UP:
140 return logFloor + lessThanBranchFree(floorPow, x);
141 case HALF_DOWN:
142 case HALF_UP:
143 case HALF_EVEN:
144 // sqrt(10) is irrational, so log10(x)-logFloor is never exactly 0.5
145 return logFloor + lessThanBranchFree(halfPowersOf10[logFloor], x);
146 default:
147 throw new AssertionError();
148 }
149 }
150
151 @GwtIncompatible("TODO")
152 static int log10Floor(long x) {
153 /*
154 * Based on Hacker's Delight Fig. 11-5, the two-table-lookup, branch-free implementation.
155 *
156 * The key idea is that based on the number of leading zeros (equivalently, floor(log2(x))),
157 * we can narrow the possible floor(log10(x)) values to two. For example, if floor(log2(x))
158 * is 6, then 64 <= x < 128, so floor(log10(x)) is either 1 or 2.
159 */
160 int y = maxLog10ForLeadingZeros[Long.numberOfLeadingZeros(x)];
161 /*
162 * y is the higher of the two possible values of floor(log10(x)). If x < 10^y, then we want the
163 * lower of the two possible values, or y - 1, otherwise, we want y.
164 */
165 return y - lessThanBranchFree(x, powersOf10[y]);
166 }
167
168 // maxLog10ForLeadingZeros[i] == floor(log10(2^(Long.SIZE - i)))
169 @VisibleForTesting static final byte[] maxLog10ForLeadingZeros = {
170 19, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12,
171 12, 12, 11, 11, 11, 10, 10, 10, 9, 9, 9, 9, 8, 8, 8, 7, 7, 7, 6, 6, 6, 6, 5, 5, 5, 4, 4, 4,
172 3, 3, 3, 3, 2, 2, 2, 1, 1, 1, 0, 0, 0 };
173
174 @GwtIncompatible("TODO")
175 @VisibleForTesting
176 static final long[] powersOf10 = {
177 1L,
178 10L,
179 100L,
180 1000L,
181 10000L,
182 100000L,
183 1000000L,
184 10000000L,
185 100000000L,
186 1000000000L,
187 10000000000L,
188 100000000000L,
189 1000000000000L,
190 10000000000000L,
191 100000000000000L,
192 1000000000000000L,
193 10000000000000000L,
194 100000000000000000L,
195 1000000000000000000L
196 };
197
198 // halfPowersOf10[i] = largest long less than 10^(i + 0.5)
199 @GwtIncompatible("TODO")
200 @VisibleForTesting
201 static final long[] halfPowersOf10 = {
202 3L,
203 31L,
204 316L,
205 3162L,
206 31622L,
207 316227L,
208 3162277L,
209 31622776L,
210 316227766L,
211 3162277660L,
212 31622776601L,
213 316227766016L,
214 3162277660168L,
215 31622776601683L,
216 316227766016837L,
217 3162277660168379L,
218 31622776601683793L,
219 316227766016837933L,
220 3162277660168379331L
221 };
222
223 /**
224 * Returns {@code b} to the {@code k}th power. Even if the result overflows, it will be equal to
225 * {@code BigInteger.valueOf(b).pow(k).longValue()}. This implementation runs in {@code O(log k)}
226 * time.
227 *
228 * @throws IllegalArgumentException if {@code k < 0}
229 */
230 @GwtIncompatible("TODO")
231 public static long pow(long b, int k) {
232 checkNonNegative("exponent", k);
233 if (-2 <= b && b <= 2) {
234 switch ((int) b) {
235 case 0:
236 return (k == 0) ? 1 : 0;
237 case 1:
238 return 1;
239 case (-1):
240 return ((k & 1) == 0) ? 1 : -1;
241 case 2:
242 return (k < Long.SIZE) ? 1L << k : 0;
243 case (-2):
244 if (k < Long.SIZE) {
245 return ((k & 1) == 0) ? 1L << k : -(1L << k);
246 } else {
247 return 0;
248 }
249 default:
250 throw new AssertionError();
251 }
252 }
253 for (long accum = 1;; k >>= 1) {
254 switch (k) {
255 case 0:
256 return accum;
257 case 1:
258 return accum * b;
259 default:
260 accum *= ((k & 1) == 0) ? 1 : b;
261 b *= b;
262 }
263 }
264 }
265
266 /**
267 * Returns the square root of {@code x}, rounded with the specified rounding mode.
268 *
269 * @throws IllegalArgumentException if {@code x < 0}
270 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and
271 * {@code sqrt(x)} is not an integer
272 */
273 @GwtIncompatible("TODO")
274 @SuppressWarnings("fallthrough")
275 public static long sqrt(long x, RoundingMode mode) {
276 checkNonNegative("x", x);
277 if (fitsInInt(x)) {
278 return IntMath.sqrt((int) x, mode);
279 }
280 /*
281 * Let k be the true value of floor(sqrt(x)), so that
282 *
283 * k * k <= x < (k + 1) * (k + 1)
284 * (double) (k * k) <= (double) x <= (double) ((k + 1) * (k + 1))
285 * since casting to double is nondecreasing.
286 * Note that the right-hand inequality is no longer strict.
287 * Math.sqrt(k * k) <= Math.sqrt(x) <= Math.sqrt((k + 1) * (k + 1))
288 * since Math.sqrt is monotonic.
289 * (long) Math.sqrt(k * k) <= (long) Math.sqrt(x) <= (long) Math.sqrt((k + 1) * (k + 1))
290 * since casting to long is monotonic
291 * k <= (long) Math.sqrt(x) <= k + 1
292 * since (long) Math.sqrt(k * k) == k, as checked exhaustively in
293 * {@link LongMathTest#testSqrtOfPerfectSquareAsDoubleIsPerfect}
294 */
295 long guess = (long) Math.sqrt(x);
296 // Note: guess is always <= FLOOR_SQRT_MAX_LONG.
297 long guessSquared = guess * guess;
298 // Note (2013-2-26): benchmarks indicate that, inscrutably enough, using if statements is
299 // faster here than using lessThanBranchFree.
300 switch (mode) {
301 case UNNECESSARY:
302 checkRoundingUnnecessary(guessSquared == x);
303 return guess;
304 case FLOOR:
305 case DOWN:
306 if (x < guessSquared) {
307 return guess - 1;
308 }
309 return guess;
310 case CEILING:
311 case UP:
312 if (x > guessSquared) {
313 return guess + 1;
314 }
315 return guess;
316 case HALF_DOWN:
317 case HALF_UP:
318 case HALF_EVEN:
319 long sqrtFloor = guess - ((x < guessSquared) ? 1 : 0);
320 long halfSquare = sqrtFloor * sqrtFloor + sqrtFloor;
321 /*
322 * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both
323 * x and halfSquare are integers, this is equivalent to testing whether or not x <=
324 * halfSquare. (We have to deal with overflow, though.)
325 *
326 * If we treat halfSquare as an unsigned long, we know that
327 * sqrtFloor^2 <= x < (sqrtFloor + 1)^2
328 * halfSquare - sqrtFloor <= x < halfSquare + sqrtFloor + 1
329 * so |x - halfSquare| <= sqrtFloor. Therefore, it's safe to treat x - halfSquare as a
330 * signed long, so lessThanBranchFree is safe for use.
331 */
332 return sqrtFloor + lessThanBranchFree(halfSquare, x);
333 default:
334 throw new AssertionError();
335 }
336 }
337
338 /**
339 * Returns the result of dividing {@code p} by {@code q}, rounding using the specified
340 * {@code RoundingMode}.
341 *
342 * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
343 * is not an integer multiple of {@code b}
344 */
345 @GwtIncompatible("TODO")
346 @SuppressWarnings("fallthrough")
347 public static long divide(long p, long q, RoundingMode mode) {
348 checkNotNull(mode);
349 long div = p / q; // throws if q == 0
350 long rem = p - q * div; // equals p % q
351
352 if (rem == 0) {
353 return div;
354 }
355
356 /*
357 * Normal Java division rounds towards 0, consistently with RoundingMode.DOWN. We just have to
358 * deal with the cases where rounding towards 0 is wrong, which typically depends on the sign of
359 * p / q.
360 *
361 * signum is 1 if p and q are both nonnegative or both negative, and -1 otherwise.
362 */
363 int signum = 1 | (int) ((p ^ q) >> (Long.SIZE - 1));
364 boolean increment;
365 switch (mode) {
366 case UNNECESSARY:
367 checkRoundingUnnecessary(rem == 0);
368 // fall through
369 case DOWN:
370 increment = false;
371 break;
372 case UP:
373 increment = true;
374 break;
375 case CEILING:
376 increment = signum > 0;
377 break;
378 case FLOOR:
379 increment = signum < 0;
380 break;
381 case HALF_EVEN:
382 case HALF_DOWN:
383 case HALF_UP:
384 long absRem = abs(rem);
385 long cmpRemToHalfDivisor = absRem - (abs(q) - absRem);
386 // subtracting two nonnegative longs can't overflow
387 // cmpRemToHalfDivisor has the same sign as compare(abs(rem), abs(q) / 2).
388 if (cmpRemToHalfDivisor == 0) { // exactly on the half mark
389 increment = (mode == HALF_UP | (mode == HALF_EVEN & (div & 1) != 0));
390 } else {
391 increment = cmpRemToHalfDivisor > 0; // closer to the UP value
392 }
393 break;
394 default:
395 throw new AssertionError();
396 }
397 return increment ? div + signum : div;
398 }
399
400 /**
401 * Returns {@code x mod m}, a non-negative value less than {@code m}.
402 * This differs from {@code x % m}, which might be negative.
403 *
404 * <p>For example:
405 *
406 * <pre> {@code
407 *
408 * mod(7, 4) == 3
409 * mod(-7, 4) == 1
410 * mod(-1, 4) == 3
411 * mod(-8, 4) == 0
412 * mod(8, 4) == 0}</pre>
413 *
414 * @throws ArithmeticException if {@code m <= 0}
415 * @see <a href="http://docs.oracle.com/javase/specs/jls/se7/html/jls-15.html#jls-15.17.3">
416 * Remainder Operator</a>
417 */
418 @GwtIncompatible("TODO")
419 public static int mod(long x, int m) {
420 // Cast is safe because the result is guaranteed in the range [0, m)
421 return (int) mod(x, (long) m);
422 }
423
424 /**
425 * Returns {@code x mod m}, a non-negative value less than {@code m}.
426 * This differs from {@code x % m}, which might be negative.
427 *
428 * <p>For example:
429 *
430 * <pre> {@code
431 *
432 * mod(7, 4) == 3
433 * mod(-7, 4) == 1
434 * mod(-1, 4) == 3
435 * mod(-8, 4) == 0
436 * mod(8, 4) == 0}</pre>
437 *
438 * @throws ArithmeticException if {@code m <= 0}
439 * @see <a href="http://docs.oracle.com/javase/specs/jls/se7/html/jls-15.html#jls-15.17.3">
440 * Remainder Operator</a>
441 */
442 @GwtIncompatible("TODO")
443 public static long mod(long x, long m) {
444 if (m <= 0) {
445 throw new ArithmeticException("Modulus must be positive");
446 }
447 long result = x % m;
448 return (result >= 0) ? result : result + m;
449 }
450
451 /**
452 * Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if
453 * {@code a == 0 && b == 0}.
454 *
455 * @throws IllegalArgumentException if {@code a < 0} or {@code b < 0}
456 */
457 public static long gcd(long a, long b) {
458 /*
459 * The reason we require both arguments to be >= 0 is because otherwise, what do you return on
460 * gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't
461 * an int.
462 */
463 checkNonNegative("a", a);
464 checkNonNegative("b", b);
465 if (a == 0) {
466 // 0 % b == 0, so b divides a, but the converse doesn't hold.
467 // BigInteger.gcd is consistent with this decision.
468 return b;
469 } else if (b == 0) {
470 return a; // similar logic
471 }
472 /*
473 * Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm.
474 * This is >60% faster than the Euclidean algorithm in benchmarks.
475 */
476 int aTwos = Long.numberOfTrailingZeros(a);
477 a >>= aTwos; // divide out all 2s
478 int bTwos = Long.numberOfTrailingZeros(b);
479 b >>= bTwos; // divide out all 2s
480 while (a != b) { // both a, b are odd
481 // The key to the binary GCD algorithm is as follows:
482 // Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b).
483 // But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two.
484
485 // We bend over backwards to avoid branching, adapting a technique from
486 // http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax
487
488 long delta = a - b; // can't overflow, since a and b are nonnegative
489
490 long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1));
491 // equivalent to Math.min(delta, 0)
492
493 a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b)
494 // a is now nonnegative and even
495
496 b += minDeltaOrZero; // sets b to min(old a, b)
497 a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b
498 }
499 return a << min(aTwos, bTwos);
500 }
501
502 /**
503 * Returns the sum of {@code a} and {@code b}, provided it does not overflow.
504 *
505 * @throws ArithmeticException if {@code a + b} overflows in signed {@code long} arithmetic
506 */
507 @GwtIncompatible("TODO")
508 public static long checkedAdd(long a, long b) {
509 long result = a + b;
510 checkNoOverflow((a ^ b) < 0 | (a ^ result) >= 0);
511 return result;
512 }
513
514 /**
515 * Returns the difference of {@code a} and {@code b}, provided it does not overflow.
516 *
517 * @throws ArithmeticException if {@code a - b} overflows in signed {@code long} arithmetic
518 */
519 @GwtIncompatible("TODO")
520 public static long checkedSubtract(long a, long b) {
521 long result = a - b;
522 checkNoOverflow((a ^ b) >= 0 | (a ^ result) >= 0);
523 return result;
524 }
525
526 /**
527 * Returns the product of {@code a} and {@code b}, provided it does not overflow.
528 *
529 * @throws ArithmeticException if {@code a * b} overflows in signed {@code long} arithmetic
530 */
531 @GwtIncompatible("TODO")
532 public static long checkedMultiply(long a, long b) {
533 // Hacker's Delight, Section 2-12
534 int leadingZeros = Long.numberOfLeadingZeros(a) + Long.numberOfLeadingZeros(~a)
535 + Long.numberOfLeadingZeros(b) + Long.numberOfLeadingZeros(~b);
536 /*
537 * If leadingZeros > Long.SIZE + 1 it's definitely fine, if it's < Long.SIZE it's definitely
538 * bad. We do the leadingZeros check to avoid the division below if at all possible.
539 *
540 * Otherwise, if b == Long.MIN_VALUE, then the only allowed values of a are 0 and 1. We take
541 * care of all a < 0 with their own check, because in particular, the case a == -1 will
542 * incorrectly pass the division check below.
543 *
544 * In all other cases, we check that either a is 0 or the result is consistent with division.
545 */
546 if (leadingZeros > Long.SIZE + 1) {
547 return a * b;
548 }
549 checkNoOverflow(leadingZeros >= Long.SIZE);
550 checkNoOverflow(a >= 0 | b != Long.MIN_VALUE);
551 long result = a * b;
552 checkNoOverflow(a == 0 || result / a == b);
553 return result;
554 }
555
556 /**
557 * Returns the {@code b} to the {@code k}th power, provided it does not overflow.
558 *
559 * @throws ArithmeticException if {@code b} to the {@code k}th power overflows in signed
560 * {@code long} arithmetic
561 */
562 @GwtIncompatible("TODO")
563 public static long checkedPow(long b, int k) {
564 checkNonNegative("exponent", k);
565 if (b >= -2 & b <= 2) {
566 switch ((int) b) {
567 case 0:
568 return (k == 0) ? 1 : 0;
569 case 1:
570 return 1;
571 case (-1):
572 return ((k & 1) == 0) ? 1 : -1;
573 case 2:
574 checkNoOverflow(k < Long.SIZE - 1);
575 return 1L << k;
576 case (-2):
577 checkNoOverflow(k < Long.SIZE);
578 return ((k & 1) == 0) ? (1L << k) : (-1L << k);
579 default:
580 throw new AssertionError();
581 }
582 }
583 long accum = 1;
584 while (true) {
585 switch (k) {
586 case 0:
587 return accum;
588 case 1:
589 return checkedMultiply(accum, b);
590 default:
591 if ((k & 1) != 0) {
592 accum = checkedMultiply(accum, b);
593 }
594 k >>= 1;
595 if (k > 0) {
596 checkNoOverflow(b <= FLOOR_SQRT_MAX_LONG);
597 b *= b;
598 }
599 }
600 }
601 }
602
603 @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L;
604
605 /**
606 * Returns {@code n!}, that is, the product of the first {@code n} positive
607 * integers, {@code 1} if {@code n == 0}, or {@link Long#MAX_VALUE} if the
608 * result does not fit in a {@code long}.
609 *
610 * @throws IllegalArgumentException if {@code n < 0}
611 */
612 @GwtIncompatible("TODO")
613 public static long factorial(int n) {
614 checkNonNegative("n", n);
615 return (n < factorials.length) ? factorials[n] : Long.MAX_VALUE;
616 }
617
618 static final long[] factorials = {
619 1L,
620 1L,
621 1L * 2,
622 1L * 2 * 3,
623 1L * 2 * 3 * 4,
624 1L * 2 * 3 * 4 * 5,
625 1L * 2 * 3 * 4 * 5 * 6,
626 1L * 2 * 3 * 4 * 5 * 6 * 7,
627 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
628 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
629 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
630 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
631 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
632 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
633 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
634 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
635 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
636 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
637 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
638 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
639 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
640 };
641
642 /**
643 * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
644 * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
645 *
646 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
647 */
648 public static long binomial(int n, int k) {
649 checkNonNegative("n", n);
650 checkNonNegative("k", k);
651 checkArgument(k <= n, "k (%s) > n (%s)", k, n);
652 if (k > (n >> 1)) {
653 k = n - k;
654 }
655 switch (k) {
656 case 0:
657 return 1;
658 case 1:
659 return n;
660 default:
661 if (n < factorials.length) {
662 return factorials[n] / (factorials[k] * factorials[n - k]);
663 } else if (k >= biggestBinomials.length || n > biggestBinomials[k]) {
664 return Long.MAX_VALUE;
665 } else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) {
666 // guaranteed not to overflow
667 long result = n--;
668 for (int i = 2; i <= k; n--, i++) {
669 result *= n;
670 result /= i;
671 }
672 return result;
673 } else {
674 int nBits = LongMath.log2(n, RoundingMode.CEILING);
675
676 long result = 1;
677 long numerator = n--;
678 long denominator = 1;
679
680 int numeratorBits = nBits;
681 // This is an upper bound on log2(numerator, ceiling).
682
683 /*
684 * We want to do this in long math for speed, but want to avoid overflow. We adapt the
685 * technique previously used by BigIntegerMath: maintain separate numerator and
686 * denominator accumulators, multiplying the fraction into result when near overflow.
687 */
688 for (int i = 2; i <= k; i++, n--) {
689 if (numeratorBits + nBits < Long.SIZE - 1) {
690 // It's definitely safe to multiply into numerator and denominator.
691 numerator *= n;
692 denominator *= i;
693 numeratorBits += nBits;
694 } else {
695 // It might not be safe to multiply into numerator and denominator,
696 // so multiply (numerator / denominator) into result.
697 result = multiplyFraction(result, numerator, denominator);
698 numerator = n;
699 denominator = i;
700 numeratorBits = nBits;
701 }
702 }
703 return multiplyFraction(result, numerator, denominator);
704 }
705 }
706 }
707
708 /**
709 * Returns (x * numerator / denominator), which is assumed to come out to an integral value.
710 */
711 static long multiplyFraction(long x, long numerator, long denominator) {
712 if (x == 1) {
713 return numerator / denominator;
714 }
715 long commonDivisor = gcd(x, denominator);
716 x /= commonDivisor;
717 denominator /= commonDivisor;
718 // We know gcd(x, denominator) = 1, and x * numerator / denominator is exact,
719 // so denominator must be a divisor of numerator.
720 return x * (numerator / denominator);
721 }
722
723 /*
724 * binomial(biggestBinomials[k], k) fits in a long, but not
725 * binomial(biggestBinomials[k] + 1, k).
726 */
727 static final int[] biggestBinomials =
728 {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 3810779, 121977, 16175, 4337, 1733,
729 887, 534, 361, 265, 206, 169, 143, 125, 111, 101, 94, 88, 83, 79, 76, 74, 72, 70, 69, 68,
730 67, 67, 66, 66, 66, 66};
731
732 /*
733 * binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl,
734 * but binomial(biggestSimpleBinomials[k] + 1, k) does.
735 */
736 @VisibleForTesting static final int[] biggestSimpleBinomials =
737 {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 2642246, 86251, 11724, 3218, 1313,
738 684, 419, 287, 214, 169, 139, 119, 105, 95, 87, 81, 76, 73, 70, 68, 66, 64, 63, 62, 62,
739 61, 61, 61};
740 // These values were generated by using checkedMultiply to see when the simple multiply/divide
741 // algorithm would lead to an overflow.
742
743 static boolean fitsInInt(long x) {
744 return (int) x == x;
745 }
746
747 /**
748 * Returns the arithmetic mean of {@code x} and {@code y}, rounded toward
749 * negative infinity. This method is resilient to overflow.
750 *
751 * @since 14.0
752 */
753 public static long mean(long x, long y) {
754 // Efficient method for computing the arithmetic mean.
755 // The alternative (x + y) / 2 fails for large values.
756 // The alternative (x + y) >>> 1 fails for negative values.
757 return (x & y) + ((x ^ y) >> 1);
758 }
759
760 private LongMath() {}
761 }